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Summary 

We describe the remapped particle-mesh method, a new mass-conserving method for solving the 
density equation which is suitable for combining with semi-Lagrangian methods for compressible flow 
applied to numerical weather prediction. In addition to the conservation property, the remapped particle- 
mesh method is computationally efficient and at least as accurate as current semi-Lagrangian methods 
based on cubic interpolation. We provide results of tests of the method in the plane, results from 
incorporating the advection method into a semi-Lagrangian method for the rotating shallow-water 
equations in planar geometry, and results from extending the method to the surface of a sphere. 
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1. Introduction 

The semi-implicit semi-Lagrangian (SISL) method, as originally introduced 
by Robert [16], has become very popular in numerical weather prediction 
(NWP). The semi-Lagrangian aspect of SISL schemes allows for a relatively accu- 
rate treatment of advection while at the same time avoiding step size restrictions 
of explicit Eulerian methods. The standard semi-Lagrangian algorithm (see, e.g., 
[19]) calculates departure points, i.e., the positions of Lagrangian particles which 
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will be advected onto the grid during the time step. The momentum and density 
equations are then solved along the trajectory of the particles. This calculation 
requires interpolation to obtain velocity and density values at the departure point. 
It has been found that cubic interpolation is both accurate and computationally 
tractable (see, e.g., [19]). 

Ideally, as well as being efficient and accurate, a density advection scheme 
should exactly preserve mass in order to be useful for, e.g., climate prediction 
or atmospheric chemistry calculations. Recent developments have involved com- 
puting the change in volume elements, defined between departure and arrival 
points, making use of a technique called cascade interpolation [14]. Several such 
methods have been suggested in recent years, including the methods of Nair et 
al [11, 12, 13] and the SLICE schemes of Zerroukat et al [23, 24, 26, 25]. 

In this paper we give a new density advection scheme, the remapped particle- 
mesh method, which is based on the particle- mesh discretisation for the density 
equation used in the Hamiltonian Particle-Mesh (HPM) method suggested by 
Gottwald, Frank & Reich [3], which itself was a combination of smoothed 
particle-hydrodynamics [7, 5] and particle-in-cell methods [6]. The particle- 
mesh method provides a very simple discretisation which conserves mass by 
construction, and may be adapted to nonplanar geometries such as the sphere 
[4]. In this paper we show that an efficient scheme can be obtained by mapping 
the particles back to the grid after each time step. Our numerical results show 
that this scheme is at least as accurate as standard semi-Lagrangian advection 
using cubic interpolation at departure points. We show how the method may be 
included in the staggered semi-Lagrangian schemes, proposed by Staniforth et 
al [20] and Reich [15], and show how to adapt it to spherical geometry. 

In section 2 we describe the particle-mesh discretisation for the density 
equation. The method is modified to form the remapped particle-mesh method in 
section 3. We discuss issues of efficient implementation in section 4. In section 6 we 
give numerical results for advection tests in planar geometry and on the sphere, 
as well as results from rotating shallow-water simulations using the remapped 
particle-mesh method in the staggered leapfrog scheme [15]. We give a summary 
of our results and discussion in section 7. 



In this section we describe the particle-mesh discretisation for the density 
equation. This discretisation forms the basis for the remapped particle-mesh 
method discussed in this paper. For simplicity, we restrict the discussion to two- 
dimensional flows. 

We begin with the continuity equation 



where p is the density and u = (u, v) T 6 M? is the fluid velocity. We write (1) in 
the Lagrangian formulation as 



2. Continuity equation and particle advection 



pt + V • (pu) = 0, 



(1) 



(2) 



Dt 
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where p(x, t) is the density at time t ^ at a fixed Eulerian position x = (x, y) T € 

R 2 , 

^(■) = (")* + (Ox u +(■)»« (4) 

is the Lagrangian time derivative, 

X(a,t) = (I(a,t),y(a,t)) T eM 2 (5) 

is a Lagrangian particle position at time t with initial position X(a, 0) = a £ M 2 , 
and po(&) = p(a, 0) is the initial density at a = (a, b) T E R 2 . 

To discretise the integral representation (3), we introduce a finite set of La- 
grangian particles X /3 (t) = (Xp(i), Yp(t)) T £ R 2 , (3 = 1, . . . , N, and a fixed Eule- 
rian grid Xkj = (xfc, yi) = (k ■ Ax, I ■ Ay) T , k, I = 0, . . . , M. Then we approximate 
the Eulerian grid density pk,i(t) pa p(x-k,l, t) by 

Pfc,j(<) : = E ^.jCM*)) dA (*fi), (6) 

where Vfc,/( x ) ^ are basis functions, which satisfy J ^^(x) &A(x) = 1. The 
initial particle positions X j a(0) = a^ are assumed to form a grid and dA(&p) is 
equal to the area of the associated grid cell. Equation (6) may be simplified to 

Pk,i(t) =5Z m / 3 ^( X /8(*))' ( 7 ) 

where 

n»/3 :=Po(a/3) dA(ap) (8) 

is the "mass" of particle /3. 

Let us now also request that the basis functions ipki satisfy the partition-of- 
unity (PolJ) property 

V> fc ,i( x ) (L4(x fc)J ) = 1, dA(x fc> ,) := AxAy, (9) 

k,l 

for all x6K 2 . This ensures that the total mass is conserved since 

Pk,l(t) dA (*k,i) = Y.Y, m P VvCMO) cb4(x M ) = mp, (10) 

k,l k,l (3 13 

which is constant. The time evolution of the particle positions X / g(t) is simply 
given by 

|*-«, (11) 

Given a time-dependent (Eulerian) velocity field u(x, t), we can discretise (8) 
and (11) in time with a simple differencing method: 

X^ +1 = X^ + Atu; +1/2 , u; +1/2 :=u(X^,t„ +1/2 ), (12) 
Pi? = E^^,/(X^ +1 ). (13) 
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In [3], this discretisation was combined with a time stepping method for the 
momentum equation to form a Hamiltonian particle-mesh method for the rotating 
shallow-water equations. The masses mp were kept constant throughout the 
simulation. In this paper, we instead combine the discretisation with a remapping 
technique so that the particles trajectories start from grid points at the beginning 
of each time step. Our remapping approach requires the assignment of new 
particle "masses" in each time step and, hence, is fundamentally different from 
semi-Lagrangian remapping strategies described, for example, in [11]. 



3. Remapped particle-mesh method 

In this section, we describe the remapped particle-mesh method for solving 
the continuity equation. The aim is to exploit the mass conservation property 
of the particle-mesh method whilst keeping an Eulerian grid data structure for 
velocity updates. To achieve this we reset the particles to an Eulerian grid point 
at the beginning of each time step, i.e., 

:= a/3 = x fc)J , = 1 + k + I ■ M. (14) 

This step requires the calculation of new particle "masses" mjg, (3 = 1, . . . , N, 
according to 

p*,j = ( 15 ) 

for given densities p^ t . This is the remapping step. We finally step the particles 
forward and calculate the new density on the Eulerian grid using equations (12)- 
(13) with nip = rnJl. Note that the Lagrangian trajectory calculation (12) can be 
replaced by any other consistent upstream approximation. Exact trajectories for 
a given time-independent velocity field u(x) will, for example, be used in the 
numerical experiments. 

The whole process is mass conserving since the PoU property (9) ensures 

that 

£ p^ 1 d^x M ) = £ £ ^ wxy- 1 ) dA(^) = £ = £ ph dA ^- 

k,l k,l (3 (3 k,l 

(16) 



4. Efficient implementation 

This density advection scheme can be made efficient since all the interpolation 
takes place on the grid; this means that the same linear system of equations, 
characterized by (15), is solved at each time step. The particle trajectories are 
uncoupled and thus may even be calculated in parallel. 

The computation of the particle masses in (15) leads to the solution of a 
sparse matrix system. We discuss this issue in detail for (area-weighted) tensor 
product cubic i?-spline basis functions, defined by 

^ /(X) := AxAy ^ cs (^A?) • ^ cs (^T ) ' (17) 
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where Y> cs (r) is the cubic B-splinc 

ipcs(r) = 

The basis functions satisfy 



| - \r\ 2 + i|r| 3 , \r\ < 1, 

|(2 -|r|) 3 , K\r\^2, (18) 

0, |r|>2. 



^^ M ( x )d^4(x M ) = l (19) 

and 

iP k ^)dA(x) = l (20) 



/ 



as required. 

A few basic manipulations reveal that (15) becomes equivalent to 

p n kjl <L4(x w ) = ^, AxAy =(l + ^fs^j (l + ^fs^j m\ x (21) 

where 

A 2 n _ ^ +li ,-2m^ + mg_ M ^Ui-H + ^-i , . 

d x m fc,/ - ^2 ' °V m k,l - ^2 ' y Zl > 

are the standard second-order central difference approximations, and we replaced 
index /3 = 1 + k + I ■ M by k, I, i.e., we write m^ l} XJ?j, etc. from now on. Eq. (21) 
implies that the particle masses can be found by solving a tridiagonal system along 
each grid line (in each direction). 

If the cubic spline Y> cs in (17) is replaced by the linear spline 

*-(»•)={ Jr |r1, t\tl: w 

then the system (15) is solved by 

m^ = AxAy^. (24) 

The resulting low-order advection scheme possesses the desirable property that 
Pki^O for all k, I implies that p^ 1 ^ for all k, I, and so that monotonicity is 
also preserved. 

On a more abstract level, conservative advection schemes can be derived 
for general (e.g. triangular) meshes with basis functions ^(x) ^ 0, which form 
a partition of unity. An appropriate quadrature formula for (3) leads then to 
a discrete approximation of type (7). This extension will be the subject of a 
forthcoming publication. 



5. Extension to the sphere 

In this section we suggest a possible implementation of the remapped particle- 
mesh method for the density equation on the sphere. The method follows the 
particle-mesh discretisation given by Frank & Reich [4], combined with a 
remapping to the grid. 
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We introduce a longitude-latitude grid with equal grid spacing AA = AO = 
tt/J. The latitude grid points are offset a half-grid length from the poles. 
Hence we obtain grid points (A&, 0\), where A& = kAX, Oi = — ^ + (I — 1/2) AO, 
k = 1, . . . , 2 J, I = 1, . . . , J, and the grid dimension is 2 J x J. 

Let ipk,i( x ) denote the (area- weighted) tensor product cubic B-spline centered 
at a grid point G R 3 with longitude-latitude coordinates (A&, 0\), i.e. 

^ (x) := di^y ^ cs (^ir) • ^ (^r) - (25) 

where (A, 0) are the spherical coordinates of a point x = (x, y, z) T G R 3 on the 
sphere, tpcs(f) is the cubic B-spline as before, and 

<L4(x fc)J ) = R 2 cos(fy) A6>AA. (26) 

We convert between Cartesian and spherical coordinates using the formulas 

x = R cos A cos 0, y = R sin A cos 0, z = Rsm0, (27) 

and 

X = tan- (|), » = (£). (28) 

At each time step we write the fluid velocity in 3D Cartesian coordinates and 
step the particles Xjj forward. We then project the particle positions onto the 
surface of the sphere as described in [4]. The Lagrangian trajectory algorithm is 
then: 

X? = x M + A< 1/2 + /1 x !J , (29) 

where p is a Lagrange multiplier chosen so that HX"^ 1 !! = R on a sphere of 
radius R. This alogrithm can be replaced by any other consistent approximation 
upstream Lagrangian trajectories. Exact trajectories are, for example, used in 
the numerical experiments. 

We compute the particle masses mfj by solving the system 

ph = T, m ijM*ij) (3°) 

for given densities p^ t . The density at time- level t n+ \ is then determined by 

p n kT=Jl m ljM x lV)- ( 31 ) 

Note that the system (30) is equivalent to 

pit dA(x fc ,,) = (l + ^fsl) (l + ml, (32) 

and can be solved efficiently as outlined in section 4. The implementation of the 
remapping method is greatly simplified by making use of the periodicity of the 
spherical coordinate system in the following sense. The periodicity is trivial in 
the longitudinal direction. For the latitude, a great circle meridian is formed by 
connecting the latitude data separated by an angular distance tt in longitude (or J 
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M 


8 


16 


32 


64 


128 


256 


512 


h 


0.549E-02 


0.254E-03 


0.143E-4 


0.872E-6 


0.541E-07 


0.337E-08 


0.211E-09 



TABLE 1. Convergence of Z2-errors as a function of Ax = 1/M for uniform advection with U = 1 of a 
sine wave on a periodic domain Q = [0, 1) with At = 0.12Ax/[/ and 20 time steps. 

grid points). See, for example, the paper by Spotz, Taylor & Swarztrauber 
[18]. It is then efficient to solve the system (32) using a direct solver. 
Conservation of mass is encoded in 

£ PI+ 1 cL^x*,,) = Pk,i dAfaj), (33) 
k,i k,i 

which holds because of the PoU property 

£^ fc ,(x)(L4(x fc) ,) = l. (34) 

k,l 

6. Numerical results 

(a) ID convergence test 

Following [26], we test the convergence rate of our method for one- 
dimensional uniform advection of a sine wave over a periodic domain fi = [0, 1). 
The initial distribution is 

po(x) = sin(27rx) (35) 

and the velocity field is u(x, t) = U = 1. The ID version of our method is used to 
solve the continuity equation 

p t = -(pu) x . (36) 

The experimental setting is equivalent to that of [26]. Table 1 displays the 
convergence of h errors as a function of resolution Ax = 1/M. Note that the 
results from Table 1 are in exact agreement with those displayed in Table I of 
[26] for the parabolic spline method (PSM) and fourth-order accuracy is observed. 

(b) 2D planar advection: Slotted- cylinder problem 

Convergence is now examined for a more realistic test case. Since we use 
higher-order interpolation the initial density profile po needs to be sufficiently 
smooth. On the other hand, relatively sharp gradients should be present to pose 
a challenge to the advection scheme. We decided to use a smoothed slotted- 
cylinder obtained by applying a modified Helmholtz operator 7i = (I — a 2 V 2 ) -1 
to the standard sharp-edged slotted cylinder [22] . The smoothing length is set to 
a = 27r/64. See panel (a) in Fig. 1. 

We compare the newly proposed scheme to the standard SL advection scheme 
based on backward trajectories and bicubic interpolation (see, e.g. [19]). To 
exclude any errors from the trajectory calculation we use a double periodic 
domain of size [0, 2ir] x [0, 2ir] and apply a constant velocity field n = 27r/3, 
v = 2tt. The time-step is At = 0.01 and the simulations are run over a period 
of T = 12 time units. Note that the initial density profile po returns to its original 
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Cubic Interpolation 
(a) initial profile 
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(b) final profile (M=256) 
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(c) difference in profiles (M=256) 
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Cubic Spline 
(d) initial profile 
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(e) final profile (M=256) 
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(f) difference in profiles (M=256) 
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Figure 1. Detailed results from the linearly advected smoothed slotted-cylindcr experiment with 
M = 256. Left panels: classic SL interpolation using backward trajectories and bicubic interpolation. 
Right panels: new advection scheme using forward trajectories and mass-conserving spline interpolation. 



position after r = 3 time units. This allows us to introduce the error 

e m HMx M )-p£f ||oo, ^ = 300, (37) 

for m = 1, 2, 3, 4. 

Simulations are performed on a spatial grid with M = 128, M = 256, and 
M = 512. Errors (37) are provided in Fig. 2. It can be seen that the newly 
proposed method is more accurate than the standard SL advection scheme and 
that the newly proposed method achieves second-order accuracy as a function 
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Figure 2. Displayed are the (oo-errors (37) for constant step-size At and varying spatial resolution 
for the lineearly advected smoothed slotted cylinder experiment. Left panel: classic SL interpolation 
using backward trajectories and bicubic interpolation. Right panel: new advection scheme using forward 
trajectories and mass-conserving spline interpolation. 



of spatial resolution (for fixed time-steps At). Detailed results from simulations 
with M = 256 can be found in Fig. 1. 

Following the discussion of [26] the reduced order can be explained by the 
fact that the Helmholtz operator H leads to an approximate \k\~ 2 spectral decay 
in the Fourier transform of the initial density pq. As also explained in [26], the 
improved convergence of our spline-based method over the traditional bicubic SL 
method is to be expected. 

We also implemented the standard rotating slotted-cylinder problem as, 
for example, defined in [10, 23]. See [23] for a detailed problem description 
and numerical reference solutions. Corresponding results for the newly proposed 
advection scheme can be found in Fig. 3. 



(c) 2D planar advection: Idealized cyclogenesis problem 
The idealized cyclogenesis problem (see, e.g., [10, 23]) consists of a circular 
vortex with a tangential velocity V(r) = vq tanh(r) /sech 2 (r), where r is the radial 
distance from the centre of the vortex (x c , y c ) and v is a constant chosen such 
that the maximum value of V(r) is unity. The analytic solution p(x, t) is 



p(x, t) = — tanh 



V ~ Vc 



cos(wi) — 



x — X, 



sin(u;i) 



(38) 



where u = V(r)/r is the angular velocity and 5 = 0.05. The experimental setting 
is that of [10, 23]. In particular, the domain of integration is SI = [0, 10] x [0, 10] 
with a 129 x 129 grid. The time step is At = 0.3125 and a total of 16 time steps 
is performed. Numerical reference solutions can be found in [23] for the standard 
bicubic and several conservative SL methods. The corresponding results for the 
newly proposed advection scheme can be found in Fig. 4. 



(d) Spherical advection: Solid body rotation 

Solid body rotation is a commonly used experiment to test an advection 
scheme over the sphere. We apply the experimental setting of [11, 12, 13, 24]. 
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Figure 3. Rotating slotted-cylindcr problem. Top panel: numerical solutation after six rotations. Bottom 
panel: error (analytic minus numerical) with contour minimum —0.5266 and contour interval 0.3803; error 
measures, as defined in [23], rmsi = 0.062595, rms 2 = 0.037329, and pdm = -0.1454E-10 %. 



The initial density is the cosine bell, 

*(M) = { Jf l 1+c ° 8 <-/*>], riR, (39) 

where R = 7ir/64, 

r = cos -1 [sin 9 + cos 9 cos(A - A c )] , (40) 

and A c = 3tt/2. The bell is advected by a time-invariant velocity field 

u = cos a cos 9 + sin a cos A sin 9, (41) 
v = — sin a sin A, (42) 

where (u, v) are the velocity components in A and 9 direction, respectively, and 
a is the angle between the axis of solid body rotation and the polar axis of the 
sphere. 

Experiments are conducted for a = 0, a = ir/2, and a = ir/2 — 0.05. Analytic 
trajectories are used and At is chosen such that 256 time steps correspond to 
a complete revolution around the globe (the radius of the sphere is set equal 
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Figure 4. Cyclogenesis problem. Top panel: numerical soluation at time t = 5. Bottom panel: error 
(analytic minus numerical) with contour minimum —0.627 and contour interval 0.418; error measures, 
as defined in [23], rmsi = 0.081439, rms 2 = 0.037703, and pdm = -0.176259E-11 %. 



a 





tt/2 


tt/2 - 0.05 


h 


0.0492 


0.0591 


0.0627 


h 


0.0336 


0.0393 


0.0397 


loo 


0.0280 


0.0367 


0.0374 



TABLE 2. Comparison ol error norms for solid body rotation with three different values of a t after 
one complete revolution using 256 time steps over a 128 X 64 grid. The meridional Courant number is 

C 9 =0.5. 



to one). Accuracy is measured as relative errors in the Zi, I2, and Zqo norms (as 
denned, for example, in [24]). Results are reported in Table 2 for a 128 x 64 grid 
(i.e., J = 64). 

Note that (32) may lead to a non-uniform distribution of particle masses 
near the polar cap regions for meridional Courant numbers Cq > 1. This can 
imply a loss of accuracy if a "heavy" extra-polar particle moves into a polar cap 
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(a) 72 time steps (b) 36 time steps (c) 18 time steps 



p 





ir/(3J) 


h 


0.0491 


0.0283 


h 


0.0468 


0.0168 


toe- 


0.0723 


0.0122 



P 





t/(3J) 


h 


2.3264 


0.0222 


h 


1.5124 


0.0137 




1.1383 


0.0151 



P 





t/(3J) 


h 


2.3217 


0.0143 


h 


1.5126 


0.0105 




1.0764 


0.0143 



TABLE 3. Comparison of error norms for solid body rotation with a = ir/2 for different values of 
the smoothing parameter (i in (43) after one complete revolution over a 128 X 64 grid (i.e., J = 64). 
Panel (a): Complete revolution using 72 time step. The meridional Courant number is Cg = 1.78. Panel 
(b): Complete revolution using 36 time step. The meridional Courant number is Cg =3.56. Panel (c): 
Complete revolution using 18 time step. The meridional Courant number is Cg = 7.12. 



region. We verified this for 72, 36 and 18, respectively, time steps per complete 
revolution (implying a meridional Courant number of Cg = 1.78, C$ = 3.56, and 
C — 9 = 7.12, respectively). It was found that the accuracy is improved by 
applying a smoothing operator along lines of constant 9 near the polar caps, 
e-g-, 

Vcos 9 3A 6 



„n+l 



n+1 



(43) 



/?<Cvr/J, J = 64. Here denotes the density approximation obtained from 

(31). The filter (43) is mass conserving and acts similarly to hyper- viscosity. The 
disadvantage of this simple filter is that p n+1 / p n under zero advection. 

Results for (3 = and (5 = 7r/192, respectively, and 72, 36 and 18 time steps, 
respectively, are reported in Table 3. It is evident that filtering by (43) improves 
the results significantly. Corresponding results for standard advection schemes 
can be found in [11] for the case of 72 time steps per complete revolution. 

(e) Spherical advection: Smooth deformational flow 

To further evaluate the accuracy of the advection scheme in spherical 
geometry, we consider the idealized vortex problem of Doswell [1]. The flow 
field is deformational and an analytic solution is available (see [9, 11] for details). 

We summarize the mathematical formulation. Let (A', 9') be a rotated 
coordinate system with the north pole at (it + 0.025, it/2.2) with respect to the 
regular spherical coordinates. We consider rotations of the (A', 6') coordinate 
system with an angular velocity u, i.e., 



It 



= to, 



dff_ 
~dt 



= 0, 



where 



cu(9') 



3\/3 sech 2 (3 cos 9') tanh(3 cos 9') 



6 cos 9' 

An analytic solution to the continuity equation (1 
provided by 

, . # „ # . „ , r3cos#' 
p(X , 9', t) = l- tanh 



(44) 



(45) 



in (A', 9') coordinates is 



sin(A' - w(0') t) 



(46) 
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t 


3 


6 


h 


0.0019 


0.0055 


h 


0.0062 


0.0172 


loo 


0.0324 


0.0792 



TABLE 4. Comparison of error norms at different times t for spherical polar vortex problem. Compu- 
tations are performed with a step size of At = 1/20 and a 128 X 64 grid. 



Simulations are performed using a 128 x 64 grid and a step size of At = 0.05. 
The filter (43) is not applied. The exact solution (evaluated over the given grid) 
and its numerical approximation at times t = 3 and t = 6 are displayed in Fig. 5. 
The relative l\, li and errors (as defined in [24]) can be found in Table 4. 
These errors are comparable to the errors reported in [11, 24] for the standard 
SL bicubic interpolation approach. 



(/) Rotating shallow-water equations in planar geometry 

To demonstrate the behavior of the new advection scheme under a time- 
dependent and compressible velocity field, we consider the shallow- water equa- 
tions (SWEs) on an /-plane [2, 17]: 



Du 

~Dt 

Dv 

~Dt 
Dfi 

~Dt 



+fv-gn x , (47) 
-fu - gfiy, (48) 
-H(ux + Vy). (49) 



Here (i = fi (x, y, t) is the fluid depth, g is the gravitational constant, and / is 
twice the (constant) angular velocity of the reference plane. 

Let H denote the maximum value of ji over the whole fluid domain. We also 
introduce the fluid depth perturbation jx = jjl — H. The perturbation satisfies the 
continuity equation 

Du _ 

— = -H (u x + v y ) (50) 

which we solve numerically using the newly proposed scheme. The overall time 
stepping procedure is given by the semi-Lagrangian Stormer-Verlet (SLSV) 
method proposed by Reich [15] with only equation (5.7) from [15] being replaced 
by the following steps: 

(i) 

^ +l/2 - £ = f ~ ^ [U X + Vy] n+1 ^ 

(ii) Solve (50) over a full time step using the newly proposed scheme with ve- 
locities (u n+1 l 2 ~ £ , v n+1 l 2 ~ e ) and initial fluid depth perturbation j] n + 1 / 2 - £ = 
^n+i/2-e _ R D enote the resulting fluid depth by /U n + 1 /2+£ = ^n+i/2+e + 
H. 
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(a) exact solution at t = 3 



(c) exact solution at t = 6 





(b) numerical solution at t = 3 




(d) numerical solution at t = 6 




Figure 5. Results of a polar vortex simulation over the sphere. The exact solution and its numerical 
approximation at time t = 3 can be found in panels (a) and (b) , respectively. Contours plotted between 
0.5 and 1.5 with contour interval 0.05. Panels (c) and (d) display the same results for time t = 6. 



(iii) 



„n+l _ „n+l/2+£ _ 



AtH 



\u x + v. 



-in+l/2+e 



The method has been implemented using the standard C-grid [2] over a 
double periodic domain with L x = L y = 3840 km (see [20] for details). The grid 
size is Ax = Ay = 60 km. The time step is At = 20 min and the value of / 
corresponds to an /-plane at 45° latitude. The reference height of the fluid is 
set to H = 9665 m. The Rossby radius of deformation is Lr ~ 3000 km. Initial 
conditions are chosen as in [20, 15] and results are displayed in an identical format 
for direct comparison. 

To assess the new discretization, results are compared to those from a two- 
time-level semi-implicit semi-Lagrangian (SISL) method with a standard bicubic 
interpolation approach to semi-Lagrangian advection (see, e.g., [8, 21]). It is 
apparent from Fig. 6 that both simulations yield similar results in terms of 
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Figure 6. Left panels: Computed time evolution, from initial time to t = 6 days, of PV over the domain 
(x, y) £ [0, 3840 km] X [0, 3840 km] using the semi-Lagrangian Stormer-Verlet (SLSV) method with time 
step At = 20 min. Contours plotted between 6.4 X 10 — 8 m~ ^-s - 1 and 2.2 X 10 — 7 m~ 1 s — 1 with contour 
interval 1.56 X 10 m s . Right panels: Differences (semi-Lagrangian Stormer-Verlet minus fully 
implicit semi-Lagrangian) at corresponding times are plotted with a 10 times smaller contour interval, 
where thin (thick) lines are positive (negative) contours. 



potential vorticity advection. Furthermore, the results displayed in Fig. 6 are 
nearly identical to those displayed in Fig. 6.1 of [15]. The implication is that 
the newly proposed advection scheme in manner very similar to the traditional 
SL interpolation scheme for this particular test problem. This result in not 
unexpected as the fluid depth remains rather smooth throughout the simulation. 
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7. Summary and outlook 

A computationally efficient and mass conserving forward trajectory semi- 
Lagrangian approach has been proposed for the solution of the continuity 
equation (1). At every time step a "mass" is assigned to each grid point which is 
then advected downstream to a (Lagrangian) position. The gridded density at the 
next time step is obtained by evaluating a bicubic spline representation with the 
advected masses as weights. The main computational cost is given by the need 
to invert tridiagonal linear systems in (21). Computationally efficient iterative or 
direct solvers are available. We also proposed an extension of the advection scheme 
to spherical geometry. A further generalization to 3D would be straightforward. 
Numerical experiments show that the new advection scheme achieves accuracy 
comparable to standard non-concerving and published conserving SL schemes. 

We note that the proposed advection scheme can be used to advect momenta 
according to 

^(pu) = -(pu)V-u. (51) 

This possibility is particularly attractive in the context of the newly proposed 
semi-Lagrangian Stormer-Verlet (SLSV) scheme [15]. 
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